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Abstract. The goal of this work is to formally abstract a Markov process evolving in 
discrete time over a general state space as a finite-state Markov chain, with the objective 
of precisely approximating its state probability distribution in time, which allows for its 
approximate, faster computation by that of the Markov chain. The approach is based 
on formal abstractions and employs an arbitrary finite partition of the state space of the 
Markov process, and the computation of average transition probabilities between partition 
sets. The abstraction technique is formal, in that it comes with guarantees on the intro¬ 
duced approximation that depend on the diameters of the partitions: as such, they can 
be tuned at will. Further in the case of Markov processes with unbounded state spaces, 
a procedure for precisely truncating the state space within a compact set is provided, to¬ 
gether with an error bound that depends on the asymptotic properties of the transition 
kernel of the original process. The overall abstraction algorithm, which practically hinges 
on piecewise constant approximations of the density functions of the Markov process, is ex¬ 
tended to higher-order function approximations: these can lead to improved error bounds 
and associated lower computational requirements. The approach is practically tested to 
compute probabilistic invariance of the Markov process under study, and is compared to 
a known alternative approach from the literature. 


1. Introduction 

Verification techniques and tools for deterministic, discrete time, finite-state systems have 
been available for many years [20j . Formal methods in the stochastic context are typically 
limited to finite-state structures, either in continuous or in discrete time UM- Stochastic 
processes evolving over continuous (uncountable) spaces are often related to undecidable 
problems (the exception being when they admit analytical solutions). It is thus of interest 
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to resort to formal approximation techniques that allow solving decidably corresponding 
problems over finite discretisations of the original models. In order to formally relate the 
computable approximate solutions to the original problems, it is of interest to come up 
with explicit bounds on the error introduced by the approximations. The use of formal 
approximations techniques over complex models can be looked at from the perspective of 
research on abstractions, which are of wide use in formal verihcation. 

Successful numerical schemes based on Markov chain approximations of general stochas¬ 
tic systems in continuous time have been introduced in the literature [2T]. However, the 
finite abstractions are only related to the original models asymptotically (at the limit, that 
is weakly), with no explicit error bounds. This approach has been applied to the approxi¬ 
mate study of probabilistic reachability or safety of stochastic hybrid models in [HIST]. An 
alternative line of work on approximations of continuous-space, discrete-time Markov pro¬ 
cesses is pursued in where the discussed approximation scheme generates a finite-state 
model. In [7] the idea of approximating by averaging is introduced, where the conditional 
expectation is used to compute the approximation, and is later extended in [6]. The weak 
point of these contributions is the fact that the approximation errors that are essential in 
assessing the quality of the approximation are not computed. 

As an alternative to qualitative approximations, in [2] a technique has been introduced 
to provide formal abstractions of discrete-time, continuous-space Markov models [Ij, with 
the objective of investigating their probabilistic invariance (safety) via model checking pro¬ 
cedures over a hnite Markov chain. In view of computational scalability, the approach has 
been improved and optimised in [iniiis], extended to a wider class of processes [EllIS], 
and practically implemented as a software tool m- These abstraction techniques hinge on 
piecewise-constant approximations of the kernels of the Markov process. Linear projection 
operators are employed in [n] to generalise these techniques via higher-order interpolations 
that provide improved error bounds on the approximation level. 

In this work we show that the approach in [21 [13] can be successfully employed to ap¬ 
proximately compute the statistics in time of a stochastic process over a continuous state 
space. We first provide a forward recursion for the approximate computation of the state 
distribution in time of the Markov process. This computation is based on a partitioning 
of the state space, and on the abstraction of the Markov process as a finite-state Markov 
chain. Further, a higher-order approximation method is presented, as a generalisation of 
the approach above, and an upper bound on the error related to the new approximation is 
formally derived. Based on the information gained from the state distribution in time, we 
show that the method can be used as an alternative to [Him in [13] to approximately com¬ 
pute probabilistic invariance (safety) for discrete-time stochastic systems over general state 
spaces. Probabilistic invariance (safety) is the dual problem to probabilistic reachability. 
Over deterministic models reachability and safety have been vastly studied in the literature, 
and computational algorithms and tools have been developed based on both forward and 
backward recursions. Similarly, for the probabilistic models under study, we compare the 
presented approach (based on forward computations) with existing approaches in the liter¬ 
ature [2l nil [T2l [T3] (which hinge on backward computations), particularly in terms of the 
introduced approximation error. 

The Markov chain abstraction applied to the forward/backward computation of proba¬ 
bilistic invariance can be generalised to other specifications expressed as non-nested PCTL 
formulae or to reward-based criteria characterised via value function recursions. Moreover, 


FORMAL ABSTRACTIONS OF CONTINUOUS-SPACE MARKOV PROCESSES 


3 


the constructed Markov chain can be shown to represent an approximate probabilistic bisim¬ 
ulation of the original process [I1E9113]. 

The article is structured as follows. Section [2] introduces the model under study and 
discusses some structural assumptions needed for the abstraction procedure. The proce¬ 
dure comprises two separate parts: Section [3] describes the truncation of the dynamics of 
the model, whereas Section [5] details the abstraction of the dynamics (approximation of the 
transition kernel) - both parts contribute to the associated approximation error. Section 
[5] considers higher-order approximation schemes and quantifies the introduced approxima¬ 
tion error. Section [6] specialises these higher-order schemes to explicit algorithms for low¬ 
dimensional models using known interpolation bases. Section [7] discusses the application of 
the procedure to the computation of probabilistic invariance, and compares it against an 
alternative approach in the literature. 

In this article we use N = {1, 2, 3,...} to denote the natural numbers, Nm = {1,2,..., m} 
and Zm = {0,1, 2, ... , m} for any m G N. 

2. Models, Preliminaries, and Goals of this work 

We consider a discrete time Markov process defined over a general state space, which is 
characterised by a pair (5, T^), where 5 is a continuous state space that we assume endowed 
with a metric and be separabl^il. We denote by {S,B{S),'P) the probability structure on 
S, with 13{S) being the Borel cr-algebrsH in S and V a probability measure to be charac¬ 
terised shortly. Tg is a stochastic kernel that assigns to each point s G 5 a probability 
measure Tg(-|s), so that for any measurable set A G B{S), P(s(l) G M|s(0) = s) = rg(^|s). 
We assume that the stochastic kernel Tg admits a density function tg, namely Tg(A|s) = 
/^tg(s|s)fis. 

Given the measurable space {S,B{S),V), we set up the product space containing 
elements s(t) = [s(0), s(l),... , s{t)], where the bold typeset is used in the sequel to indicate 
vector quantities. Suppose that the initial state of the Markov process is distributed 
according to the density function vro : 5 —?■ Then the multi-variate density function 

to('So)G('Si|so)G('S 2 |si) .. .ts{st\st-i) is a probability measure P on the product space 5*+^. 
On the other hand the state distribution of at time t G N is characterised by a density 
function vr^ : 5 —?■ R-*^, which fully describes the statistics of the process at time t and is in 
particular such that, for all M G .8(5), 

P(s(t) e A) = [ Trtis)ds, 

J A 

where the symbol P is used to indicate the probability associated to events over the product 
space 5*^^ (note that the event s{t) G ^ is equivalent to s(t) G 5* x A on their corresponding 
probability spaces). 

The state density functions 7rt(-) can be characterised recursively, as follows: 

7rt+i(s) = j ts{s\s)Trt{s)ds Vs G 5. (2.1) 

metric space S is called separable if it has a countable dense subset. 

^The Borel cr-algebra B{S) is the smallest cr-algebra in S that contains all open subsets of S. For a 
separable metric space S, B{S) equals the n-algebra generated by the open (or closed) balls of S. 
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In practice the forward recursion in (12. ip rarely yields a closed form for 7rt(-). A special 
instance where this is the case is represented by a linear dynamical system perturbed by a 
Gaussian process noise: due to the closure property of Gaussian distributions over addition 
and multiplication by a constant, it is possible to explicitly write recursive formulae for the 
mean and the variance of the distribution, and thus express in a closed form the distribution 
in time of the solution process. In more general cases, it is necessary to numerically (hence, 
approximately) compute this density function in time. 

This article provides a numerical approximation of the density function of in time as 
the probability mass function (pmf) of a finite-state Markov chain ^p. The Markov chain 
^p is obtained as an abstraction of the concrete Markov process The abstraction is 
associated with a guaranteed and tunable error bound, and algorithmically it leverages a 
state-space partitioning procedure. The procedure is comprised of two steps: 

(1) since the state space S is generally unbounded, it is first properly truncated; 

(2) subsequently, a partition of the truncated dynamics is introduced. 

Section [3] discusses the error generated by the state-space truncation, whereas Section 0] 
describes the construction of the Markov chain by state-space partitioning. The discussed 
Markov chain abstraction is based on a piecewise-constant approximation of the density 
functions. In order to improve the efficiency and the precision of the approximation, we 
generalise the abstraction method in Section [5] utilizing higher-order approximations of the 
density functions. We employ the following example throughout the article as a running 
case study. 

Example 2.1. Consider the one-dimensional stochastic dynamical system 

s(t -|- 1) = as{t) -1-6-1- aw{t), 

where the parameters a,iT > 0, whereas 6 G M, and w{-) is a process comprised of indepen¬ 
dent, identically distributed random variables with a standard normal distribution. The 
initial state of the process is selected uniformly within the bounded interval [/3o,7o] C M. 
The solution of the model is a Markov process, evolving over the state space 5 = M, and 
fully characterised by the conditional density function 

ts(s|s) = (j)a{s — as — b), where (pa{u) = —. □ 

crv2vr 

We raise the following assumptions in order to relate the state density function of 
to the probability mass function of ^p. More precisely, these assumptions are employed 
for the computation of the approximation error, but the abstraction approach proposed in 
this article can be applied without raising them. 

Assumption 2.2. For given sets F C and Aq C S, there exist positive constants e and 
Eq, such that fs(s|s) and vro(s) satisfy the following conditions: 

is(s|s) < e, V(s,s) G S‘^\F, and vro(s) < eq, Vs G 5\Ao. (2.2) 

Assumption 2.3. The density functions vro(s) and ts(s|s) are (globally) Lipschitz continu¬ 
ous, namely there exist finite constants Ao,Af, such that the following Lipschitz continuity 
conditions hold: 

|7ro(s) - 7ro(s')| < Ao||s - s'||. Vs, s' G Aq, 

|4(s|s) — ts(s'|s)| < Af||s — s'll. Vs,s,s' G S. 


(2.3) 

(2.4) 
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Moreover, there exists a finite constant such that 

Mf = sup Iy ts{s\s)ds 

The Lipschitz constants Aq, Xf are practically computed by taking partial derivatives 
of the density functions 7 ro(-), t 5 (-|s) and maximising their norm. The sets Aq and F will 
be used to truncate the support of the density functions 7 ro(-) and t 0 (-|s), respectively. 
Assumption 12.21 enables the precise study of the behaviour of density functions 7rt{-) over the 
truncated state space. Furthermore, the Lipschitz continuity conditions in Assumption 12.31 
are essential to derive error bounds related to the abstraction of the Markov process over 
the truncated state space. In order to compute these error bounds, we assign the infinity 
norm to the space of bounded measurable functions over the state space S, namely 

||/(■s)||cx> = sup|/(s)|, V/ E ]B(5) = {s' : 5 —>■ M, g bounded and measurable}. 

sGS 

In the sequel the function lyi(-) denotes the indicator function of a set ACS, namely 
Ia('S) = I, if s E A; else = 0. 

Example 12.11 (Continued). Select the interval Aq = [/So) 7 o] and define the set T by 
the linear inequality 

r = {(s, s) E M^l |s — as — 6 | < aa}. 

The initial density function ttq of the process can be represented by the function 

7ro(s) = 1 [^q_.^j,](s)/(7o -/3o)- 

Then Assumption 12.21 is valid with constants e = 4>i{a)/cr and eo = 0. The constant in 
Assumption 12.31 is equal to 1/a. Lipschitz continuity, as per (12.31) and (12.411 . holds for the 
constants Aq = 0 and Aj = 1/ (cj^\/27re )■ □ 


e5 


(2.5) 


3. State-Space Truncation Procedure 

We limit the support of the density functions ttq , t^ to the sets Aq , F respectively, and 
recursively compute support sets At, as in (13.2p . that are associated to the density functions 
TTt. Then we employ the quantities e, eo in Assumption [2]2] to compute bounds et, as in (|3.1I) . 
on the error incurring in disregarding the value of the density functions vrt outside the sets 
At. Finally we truncate the original, unbounded state space to the set T = U^QAt. 

As intuitive, the error related to the spatial truncation depends on the behaviour of 
the conditional density function tg over the eliminated regions of the state space. Suppose 
that sets F,Aq are selected such that Assumption 12.21 is satished with constants €,Sq: then 
Theorem 13.31 provides an upper bound on the error obtained from manipulating the density 
functions in time 7rt(-) exclusively over the truncated regions of the state space. 

Theorem 3.1. Under Assumption\2/M the functions nt satisfy the bound 

0 < 7Tt{s) < et, Vs E S\At, 
where the quantities {et, t E are defined recursively by 

et+i = e -|- M ^ et , (3-1) 

whereas the support sets {At, t E Ztvj are computed as 

At+i = n,- (r n (At x 5)), 


(3.2) 
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where denotes the projection map along the second set of coordinate^. 


Remark 3.2. Notice that if the shape of the sets F and Aq is computationally manageable 
{e.g., if these sets are polytopes), then it is possible to precisely implement the computation 
of the recursion in ()3.2p by available software tools, such as the MPT toolbox [22]. Further, 
if for some to, D A^^, then for all t > to, At+i D A*. Similarly, we have that 

• if for some to, A^g+i C A^p, then for all t > to, A^+i C A^. 

• if for some to, A^g+i = At^, then for all t > to, A^ = A^g. 

In order to clarify the role of F in the computation of At, we emphasize that A^+i = 
UsgAt“(s), where S depends only on F and is dehned by the set-valued map 

S:S^2^, S{s) = {s eS\{s,s) e F}. 

Figured) provides a visual illustration of the recursion in (13.2p . □ 


Let us introduce the quantity K{t,Mf), which plays a role in the solution of ()3.ip and 
will be frequently used shortly: 


K{t, Mf) 


l-Mf 
l-Mf ’ 

L 


Mf / 1 
Mf = 1 . 


(3.3) 


The following theorem provides a truncation procedure, valid over a finite time horizon 
Zat = {0,1,... , N}, which reduces the state space S to the set T = IJt^o theorem 

also formally quantihes the associated truncation error. 


Theorem 3.3. Suppose that the state space of the proeess has been truneated to the set 
T = introduee the following recursion to compute funetions p-t : S ^ 

as an approximation of the density funetions nt: 

Ht+i{s) = lr{s) j ts{s\s)nt{s)ds, nois) = 1aq{s)tto{s), Vs G 5. (3.4) 

Then the introduced approximation error is \\T^t — Tt\\oo If for all t ^'Em- 

Recapitulating, Theorem l3.3l leads to the following procedure to approximate the density 
functions of .^#3 over an unbounded state space S : 

(1) truncate tto in such a way that p,o has a bounded support Aq; 

(2) truncate the conditional density function ts(-|s) over a bounded set for all s G 5, then 
quantify T C 5^ as the support of the truncated density function; 

(3) leverage the recursion in (|3.2I) to compute the support sets Ap, 

(4) use the recursion in ()3.4I) to compute the approximate density functions p,t over the set 
T = U^^gAf. Note that the recursion in (13.41) is effectively computed over the set T, 
since /xt(s) = 0 for all s G 5\T. 

Note that we could as well handle the support of //*(•) over the time-varying sets At, by 
adapting the recursion in (|3.4I) with lAt+i instead of It- However, while employing the 
(larger) set T may lead to an increase in memory requirements at each stage, it will consid¬ 
erably simplify the computations of the state-space partitioning and of the Markov chain 
abstraction: indeed, employing time-varying sets At would render the partitioning proce¬ 
dure also time-dependent, and the obtained Markov chain would be time-inhomogeneous. 
We therefore opt to work directly with set T in order to avoid these difficulties. 


^Recall that both F and A x S are defined over = S x S. 
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Figure 1: Graphical representation of the recursion in (13.2h for sets A^. 

Example 12.11 (Continued). We easily obtain a closed form for the sets = [/3t,7t], via 

I3t+i = a(3t + b- aa, jt+i = a'yt + b + aa. 

Set T is the union of intervals [(5t , 7t] . The error of the state-space truncation over T is 

Ikt -/Utiloo < et = Mf =-. 

a a 

□ 

The recursion in (j3.4l) indicates that the support of functions ixt are always contained 
in the set T (namely, they are equal to zero over the complement of T). We are thus 
only interested in computing these functions over the set T, which allows simplifying the 
recursion in (j3.4l) as follows: 

Ht+i{s) = [ ts{s\s)ntis)ds, Vs E T, tEZAT-i. (3.5) 

Jr 

4. Markov Chain abstraction via State-Space Partitioning 

In this section we assume that sets F, Aq have been properly selected so that T = U^gA^ 
is bounded. In order to formally abstract process as a finite Markov chain and to 
approximate its state density functions, we select a finite partition of the bounded set T as 
T = where sets Ai have non-trivial measure. We then complete the partition over 

the whole state space S = by additionally including set An+i = 5\T. This results in 

a finite Markov chain with n-|-l discrete abstract states in the set = {1, 2,... , n, n-|- 
1}, and characterised by the transition probability matrix P = {Pij\ E where 

the probability of jumping from any pair of states i to j (Pij) is computed as 

~ c{Ai) fAj fAi is(sls)dsds, Vi E N^, 

for all j E Nn-i-i, and where 6 (^+ 1 )j is the Kronecker delta function (the abstract state n -|- 1 
of ^p is absorbing), and £(•) denotes the Lebesgue measure of a set (i.e., its volume). The 
quantities in (14.11) are well-defined since the set T is bounded and the measures C{Ai),i E 
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Nn, are finite and non-trivial. Notice that matrix P for the Markov chain is stochastic, 
namely 



1 

C{Ai) 



ts{s\s)dsds 


1 

C{Ai) 



ds = 1. 


The initial distribution of is the pmf po = \po{l),po{2),... ,po{n + l)], and it is obtained 
from TTo as po{i) = 7ro(s)ds, Vi G Then the pmf associated to the state distribution 

of ^p at time t can be computed as pt = po-P*. 

It is intuitive that the discrete pmf pt of the Markov chain ^p approximates the 
continuous density function vrt of the Markov process la the rest of the section we show 
how to formalise this relationship; pt is used to construct an approximate function, denoted 
by ipt, of the density Tit- Theorem 14.21 shows that ipt is a piecewise constant approximation 
(with values in its codomain that are the entries of the pmf pt, normalised by the Lebesgue 
measure of the associated partition set) of the original density function tt^. Moreover, under 
the continuity assumption in (12.4h (ref. Lemma Fd.ljl we can establish the Lipschitz continuity 
of TTf , which enables the quantification (in Theorem 14.2h of the error related to its piecewise 
constant approximation 


Lemma 4.1. Suppose that the inequality in (I23D holds. Then the state density funetions 
7rt{-) are globally Lipschitz continuous with constant Af for all t G N.' 

|7rt(s) - 7rt(s')l < Af||s - s'||. Vs, s' G S. 

Theorem 4.2. Under Assumvtions \2.2\ and \2.3\. the functions TTt{-) can be approximated by 
piecewise constant functions ipti'), defined as 

= VtGN, (4.2) 

i=i 

where 1_b(') is the indicator function of a set B G S. The approximation error is upper- 
bounded by the quantity 

WT^t - iptWoo < £t + Et, Vt G N, (4.3) 

where can be recursively computed as 

Et+i = MfEt + Af(5, Eq = Xo6, (4-4) 

and 6 is an upper bound on the diameters of the partition sets {Ai}'l^i, namely S = 
sup{||s — s'll, s,s'GAlj) i G N„}. 

Note that the functions ift are dehned over the whole state space S, but (14.2h implies 
that they are equal to zero outside the set T. 


Corollary 4.3. The recursion in (14.4p admits the explicit solution 

Et= [A(t,Mf)Af + M*Ao] 5, 

where K{t,Mf) is introduced in (13.31) . 
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Underlying Theorem 14.21 is the fact that V’t(') general sub-stochastic density 

functions: 


n / .\ n 

= = 1 + 1 ) < 1 . 
i=i *1 i=i 

This is clearly due to the fact that we are operating on the dynamics of truncated 
over the set T. It is thus intuitive that the approximation procedure and the derived error 
bounds are also valid for the case of sub-stochastic density functions [3], namely 

/ ts{s\s)ds <1, Vs G 5, / 7 ro(s)(is < 1, 

Js Js 

the only difference being that the obtained Markov chain is as well sub-stochastic. 

Further, whenever the Lipschitz continuity requirement on the initial density function, 
as per (j2.3p in Assumption 12.31 does not hold, (for instance, this is the case when the initial 
state of the process is deterministic) we can relax this continuity assumption on the initial 
distribution of the process by starting the discrete computation from the time step t = 1. 
In this case we define the pmf pi = [pi(l),pi(2),... ,pi{n + 1)], where 


Pi{i) = / / ts{s\s)7ro{s)dsds, Vf G N„+i, (4.5) 

■JAi Js 

and derive pt = piP*“^ for all t G N. Theorem 14.21 follows along similar lines, except for 
equation ( 031 ), where the initial error is set to Eq = 0 and the time-dependent terms Et 
can be derived as Et = K{t, 

It is important to emphasise the practical computability of the derived errors, and 
the fact that they can be tuned by selecting a finer partition of set T that relates to a 
smaller global parameter 6 . Further, in order to attain abstractions that are practically 
useful, it imperative to seek improvements on the derived error bounds: in particular, 
the approximation errors can be computed locally (under corresponding local Lipschitz 
continuity assumptions), following the procedures discussed in |13j . 


Example 12.11 (Continued). The error of the Markov chain abstraction can be expressed 
as 


Ft - 


< tv(t, Mf) 


+ 




1 

Mf = -. 
a 


(4.6) 


The error can be tuned in two distinct ways: 

( 1 ) by selecting larger values for a, which on the one hand leads to a less narrow truncation, 
but on the other requires the partition of a larger interval; 

( 2 ) by reducing partitions diameter 6 , which of course results in a larger cardinality of the 
partition sets. 


Let us select values 6 = 0, /3o = 0, 70 = 1, <t = 0.1, and time horizon N = 5. For a = 1.2 we 
need to partition the interval T = [—0.75a, 2.49 -|- 0.75a], which results in the error jj-TTf — 
Ft Iloo < 86.8(5-|-35.9(/>i (a) for all t G Z^v- For a = 0.8 we need to partition the smaller interval 
T = [—0.34a, 0.33 -|- 0.34a], which results in the error \\TTt — iptWoa < 198.6(5 -|- 82.14>i{a) for 
all t G Zfv- Notice that in the case of a = 1.2, we partition a larger interval and obtain 
a smaller error, while for a = 0.8 we partition a smaller interval with correspondingly a 
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Figure 2: Piecewise constant approximation ipti') of the state density function 7rt{-) (derived 
analytically), for parameters a = 1.2 (left) and a = 0.8 (right). 


larger error. It is obvious that the parameters 5, a can be chosen properly to ensure that a 
certain error precision is met. This simple model admits a solution in closed form, and its 
state density functions can be obtained as the convolution of a uniform distribution (the 
contribution of initial state) and a zero-mean Gaussian distribution with time-dependent 
variance (the contributions of the process noise). This leads to the plots in Figure [2l which 
display the original and the approximated state density functions for the set of parameters 
a = 2.4,5 = 0.05. □ 


5. Higher-Order Approximation Schemes 

In the previous section we have shown that a Markov chain abstraction can be employed to 
formally approximate the density function of a Markov process in time. This abstraction 
is interpreted as a piecewise-constant approximation of the density function of the Markov 
model. In this section we argue that this procedure can be extended to approximation 
techniques based on higher-order interpolations. 

With focus on the truncated region of the state space, let us denote with B(T) the 
space of bounded and measurable functions / ; T —)• M, equipped with the infinity norm 
ll/lloo = sup{|/(s)|, s E T}, for all / E ]B(T). The linear operator TZj;, defined over ]B(T) 
by 

nrifm = [ t,{s\s)f{s)ds, VseT, V/eB(T), (5.1) 

Jr 

characterises the solution of the recursion in (13.51) as fitis) = Tly{fj,o){s), for any t E N^r. 
While in Section[4]we have proposed approximations of functions //*(•) by piecewise-constant 
functions V’t(') with an explicit quantification of the associated error, in this section we are 
interested in considering approximations via higher-order interpolations. 
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5.1. Quantification of the Error of a Projection Over a Function Space. Consider 
a set of basis functions $ = {(^i(s), (/> 2 (s),..., h £ N, the function space ib = 

span^ generated by this set as a subset of ]B(T), and a linear operator fix : B(T) —)• 'b, 
which projects any function / E B(T) onto the function space ib. Theorem 15.11 provides 
a theoretical result for approximating the solution of (|3.5p : the following section provides 
details on turning this result into a useful tool for approximations. 


Theorem 5.1. Assume 

that a linear projection operator IIx : B(T) 

— )■ T satisfies the 

inequality 

\\Ur{tff\s))-tff\s)\\^<£\ VseT, 

(5.2) 

and that there exists a finite constant , such that 



j |nx(ts(s|s))|ds < Mj’’, VseT. 

(5.3) 

Define the functions 

) as approximations of pt{-) (cf. (15. ip ). by 



i’t = (nx'E.x)*(Aio)) t £ '^N- 

(5.4) 

Then it holds that 

\\pt - iPtWoo < , teNN, 

(5.5) 

where the error e\ satisfies the difference equation 



o 

II 

+ 

II 

bq 



Corollary 5.2. Under the assumptions raised in (I5.2l) - (l5.3p . the error e\ can he alterna¬ 
tively expressed explicitly as e\ = £^K{t,M^). 

The error e\ formulated in Theorem 15.II is comparable with the quantity Et computed 
in Theorem 14.21 Both , Et represent bounds on the approximation error introduced by 
the density function obtained after state space truncation. The difference is in the 
initialisation of the corresponding recursions, where we have Eg = 0 because i/’o = but 
Eq = Xq5 since i/jq is a piecewise constant approximation of po in (14.21) . As we mentioned 
before, the quantities in (|4.5I) can be alternatively employed as the starting values of the 
computation to relax the continuity assumption on ttq, which results in an initial error 
Eq = 0, thus providing a complete similarity between Et and e\. 

5.2. Construction of the Projection Operator. In the ensuing sections we focus, for 
the sake of simplicity, on a Enclidean domain, namely T C 5 = where d denotes a finite 
dimension. We discuss a general form for an interpolation operator related to the discussed 
projection operation. Let (j)j -.V ^ E N/i, be independent fnnctions defined over 

a generic set T>. The interpolation operator np is defined as a projection map onto the 
function space T = span{4 )i{s), (f) 2 {s ),..., 4 >h{s)}, which projects any function / : P —)• M to 
a unique function o:j4’j^ using a finite set of data {{sj, f{sj))\sj £'D,j E N/i} 

and such that Ilx>{f){sj) = f {sj). The projection coefficients aj,j E N/^, satisfy the linear 
equation f = Qa, where f = [/(sj)]jgN;i and a = are L-dimensional column vectors, 

and Q = [4>j{si)]ij is the associated {h x /i)-dimensional interpolation matrix. 

Let us now shift the focus to the recursion in ()3.5p discussed in the previous section and 
tailor the operators above accordingly. Let us select a partition {.4.*, i E N^} for the set T, 
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with finite cardinality n. Selecting a basis j G for each partition, let us introduce 
the interpolation operators n_ 4 . for the projection over each partition set Ai, which is done 
as described above by replacing the domain T> with Ai. Finally, let us introduce the (global) 
linear operator IIx, acting on a function / : T —)• M by 

n 

nT(/) = ^lAnA(/lA), (5.6) 

i=l 

where /|a represents the restriction of the domain of function / to the partition set Ai- 


5.3. Approximation Algorithm. An advantage of the interpolation operator in (j5.6|] is 
that nT(/) is fully characterised by the interpolation coefficients aij, since 

n h 

= EE 

i=l j=l 

The set of interpolation coefficients Oij is computable by matrix multiplication based on 
the data set {/(sij), i G N„, j G N/i}. More precisely, we have 

with the interpolation matrices Qi = These matrices depend solely on the 

interpolation points Sij and on the basis functions (j)ij evaluated at these points and can be 
computed off-line (see step 0] in Algorithm [H to be discussed shortly). Moreover, values of 
the function / only at the interpolation points Sij are sufficient for the computation of a^j. 

Let us now focus on the recursion in (15.4p . namely V’t+i = i^^t)i given the initial¬ 

isation ?/)g = hq, for the approximate computation of the value functions. This recursion 
indicates that the approximate functions V't > ^ G Nat, belong to the image of the operator 
fix, and as such can be expressed as 

n h 

i=i j=i 

where ajj denote the interpolation coefficients referring to function (at step t). This 
suggests that we need to store and update the coefficients ajj for each iteration in (j5.4l) . 

Writing the recursion in the form V'i+i = nx(7^x(V’t)) indicates that the function 
is in the range of the projection Fix. Therefore, it is sufficient to evaluate the function 
TZr{'4’t) over the interpolation points in order to compute the coefficients c/ij^- In the 
following expressions, the pair i, u indicates the indices of related partition sets, namely 
Ai,Au, whereas the pair of indices j,v show the ordering positions within partition sets. 
For an arbitrary interpolation point Suv we have: 


~ n h „ 

TZr(suv) = / ts{suv\s)i)t{s)ds = / ts{sst\s)c/)ij{s)ds. 

Jt JAi 


Introducing the following quantities 


puv 

^ij 



{,^uv I {s)ds^ 
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we can succinctly express 

n h 

1=1 j=i 

Algorithm [1] provides a general procedure for the discrete computation of the interpolation 
coefficients and of the approximate value functions. 


Algorithm 1 Approximate computation of the functions 
Require: Density function ts(s|s), set T 

1: Select a finite n-dimensional partition of the set T = (Ai are non-overlapping) 

2: For each A*, select interpolation basis functions (pij and points Sij G Ai, where j G 
3: Compute t^(suv\s)c/>ij(s)ds, where i,u G Nn and j,v G Nh 

4: Compute a matrix representation for the operators n_ 4 ., namely Qi = {(t>ij{siv)]v,j 
5: Set t = 1 and /3/j = for all i,j 

6: if t < A then 

7: Compute interpolation coefficients aC based on equation 

given /3b and matrices Qi in step [4] 

8: Compute values /3*+^ as /3‘+^ = Ylj ^ 

9: t = t + l 

10 : end if 

Ensure: Approximate functions V’t = Yli t £ Ntv 


It is possible to simplify Algorithm [T] when the interpolation matrices Qi are nonsingular. 
Let us transform the basis {4>ij, j G N/j} to its equivalent basis using matrix {QJ) ■ The 
interpolation matrices corresponding to the new basis will be the identity matrix. In other 
words, the new basis functions admit Qi = Ih and thus step 0] can be skipped, and that the 
main update (steps [7] and [8|) can be simplified as follows: 

n h 

= E E N;,. 

i=l j=l 

In Algorithm [H the interpolation points Sij are in general pair-wise distinct. By ex¬ 
tending the domain of interpolation Ai to its closure Ai, it is legitimate to use boundary 
points as interpolation points, which can lead to a reduction of the number of integrations 
required in Algorithm [TJ In the ensuing sections, we will exploit this feature by specifically 
selecting equally spaced interpolation points. 

6. Special Forms of the Projection Operator 

In this section we leverage known interpolation theorems for the construction of the projec¬ 
tion operator IIx: this should both yield useful schemes for a number of standard models, 
and further help with the understanding of the details discussed in the previous section. 
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6.1. Piecewise Constant Approximations. We focus on the special case of the approx¬ 
imation of a function by a piecewise constant one, which has inspired Section HI Let us 
select the basis functions = 1 for all i G N„,, j G Ni - the cardinality of these sets 

of basis functions is simply equal to h = 1 (we eliminate the corresponding indices when 
appropriate). In this case the matrix operators Qi, i G {cf. step [H in Algorithm [T]) 
correspond to the identity matrix, and the projection operator Xl-r becomes 

n 

nT(/) = 5]/(5*)lA, V/gB(T), (6.1) 

i=l 

where the quantities (cf. step [3] in Algorithm [1]) form a square matrix (see step [3] 
in Algorithm [2]) . The procedure is detailed in Algorithm [2l while the associated error is 
formally quantified in Theorem 16.11 


Algorithm 2 Piecewise constant computation of the functions 
Require: Density function ts(s|s), set T 

1: Select a finite n-dimensional partition of the set T = {Ai are non-overlapping) 

2: For each Ai, select one representative point s* G Ai 

3: Compute matrix P = [P{i,j)\i,j with entries P{i,j) = ts{sj\s)ds, where i,j G 
4: Set t = l and ai{i) = Jyts{si\s)fj:o{s)ds, for all i 
5: if t < N then 

6: Compute the row vector = [at+i{i)\i based on at+i = cxtP 

7: t = t + l 

8 : end if 

Ensure: Approximate functions i/’t = t G Nat 


Theorem 6.1. Suppose the density function ts(js) satisfies the Lipschitz continuity as¬ 
sumption with constant Ap Then the projection operator m satisfies the inequality 

||nT(t3(-|s))-ta(-|s)||oo < Afd, VsGT, 

where 6 = maxj d* is the partition diameter ofU^^^Ai = T, with 6i = sup{||s — s'|| : s,s' G 
Ai}. Theorem \5.1\ ensures that the approximation error of Algorithmic is upper bounded by 
the quantity 

WlJ-t - f^tWoo < = XfSK{t,M^), t G Nat, 

with the constant Mf defined in Assumption\2fC 

Notice that the error of Theorem 16.11 reduces to Et in Theorem 14.21 when employing 
the quantities in (|4.5p to relax the continuity assumption on tto. 

Let us compare Algorithms [T] and [2] in terms of their computational complexity. Algo¬ 
rithm [1] requires nh{nh P 1) integrations in the marginalisation steps ([3] and [5]), whereas 
n(n -|- 1) integrations are required in Algorithm [2j Furthermore, steps H] and [3 in Algo¬ 
rithm [1] can be skipped by using proper equivalent basis functions, whereas these steps are 
not needed at all in Algorithm [2j As a bottom line, higher interpolation orders increase 
the computational complexity of the approximation procedure, however this can as well 
lead to a lower global approximation error. From a different perspective, since the global 
approximation error depends on the local partitioning sets (their diameter and the local 
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continuity of the density function), for a given error higher interpolation procedures may 
require less partitions sets. 

As a hnal note, comparing the transition probabilities of (14.ip with quantities P{i,j) in 
step [3] of Algorithm [2] reveals that the Markov chain abstraction presented in Section d] is a 
special case of Algorithm [2j More precisely, the mean value theorem for integration ensures 
the existence of representative points Sj such that P{i,j) of Algorithm [2] is equal to Pij in 

gU). 


6.2. Higher-order Approximations for One-Dimensional Systems. We study higher- 
order interpolations over the real axis, where the partition sets Ai are real-valued intervals. 
We use this simple setting to quantify the error related to the approximate computation of 
the functions fit- We select equally spaced points as the interpolation points and employ 
polynomial basis functions within each interval. 

Consider a one dimensional Markov process, 5 = M, with a partitioning of T = 
which is such that Ai = [ai,bi]- Define the interpolation operator IIx of (15.6|) over the 
polynomial basis functions i G N„, j € N/^,/i > 2, (or their equivalent 

Lagrange polynomials [24]) using equally spaced interpolation points Sij G Ai, 

Sij = ai + {j - , j G Nft. 

The following result can be adapted from |24j . 


Theorem 6.2. Assume that the density function ts(-|s) is h-times differentiable and define 
the constant 


Mh 


max 

s,seT 


d^jsjs) 

ds^ 


The interpolation operator fix, constructed with polynomial basis functions and equally 
spaced interpolation points, satisfies the inequality 

IFr («,(F)) - (.(Fjlloo < ^ € T, 

where 6 = maxj(ij, with 6i = bi — ai, i G and where h is the cardinality of the set of 
basis functions. 


Theorem 16.21 provides the necessary ingredients for Theorem l5.ll leading to the quantifi¬ 
cation of the approximation error: employing Algorithm [T] with equally spaced points and 
polynomial basis functions of degree less than h, the approximation error is upper bounded 
by the quantity 

ll/it - V’t Iloo < E^t = £''K{t,Mf), t G Nn, 

with the constant defined in (15.3p and computed for this particular choice of basis 
functions and points. 

It is worth highlighting that, unlike the piecewise constant case of Section d] with higher- 
order approximation approaches the global error is a nonlinear function of the partition size 
6 , namely it depends on a power of the partition size contingent on the order of the selected 
interpolation operator. As such, its convergence speed, as 5 is decreased, increases over that 
of the piecewise constant case. 








16 


S. ESMAEIL ZADEH SOUDJANI AND A. ABATE 


a = 1.2 



a = 0.8 



Figure 3: Comparison of the first-order (affine) approximation versus the Markov 

chain abstraction (zero-order approximation, constant) '0Af(') of the state density 
function (derived analytically), for = 5 and parameters a = 1.2 (left) 

and a = 0.8 (right). 


Example 12.11 (Continued). We partition the set T = for the one dimensional 

system of Example 12.11 with the intervals Ai = [oj, bi]. We select interpolation points 
{oj, Uj+i} with polynomial basis functions {1, s}, leading to piecewise affine approximations 
(namely, first-order interpolation with h = 2) of the density functions vri(-). This set of basis 
functions can be equivalently transformed to 



to obtain Qi = II 2 - The constant has the same value as = Ija and the quantity Ad 2 
in Theorem 16.21 is AI 2 = The error related to this first-order approximation can 

be upper bounded as 


TTt - V’t Iloo < K(t,Mf) 




<(>i(«) 

CJ 


( 6 . 2 ) 


Notice that the first part of the error in ()6.2p . which specifically relates to the first-order 
approximation, is proportional to 5“^ - this improves the error bound computed in (14.61) . 
Algorithm [T] is implemented for this linear system with the aforementioned parameters 
6 = 0,/3o = 0,7o = 1,(T = 0.1, and the time horizon N = 5. The errors corresponding 
to the values a = 1.2 and a = 0.8 are analytically upper-bounded, for any t G Nat, as 
IKt “ V’t Iloo — 179(5^ -|- 35.9(/ii(a) and as Hvr^ — V't ||oo < 409.35^ -|- 82.1(/)i(a), respectively. 
The plots in Figure [3] display the first- and zero-order approximations of the density function 
7r7v(-) and compare it with the analytic solution for two different values a = 1.2 (left) and 
a = 0.8 (right). The partition size n = 25 and parameter a. = 2.4 have been selected in 
order to illustrate the differences of the two approximation methods in Figure [3l but may 
be increased at will in order to decrease the related error bound to match a desired value. 


6.3. Bilinear Interpolation for Two-Dimensional Systems. We directly tailor the 
results of Section [5]to a general two-dimensional Markov process, where s = (si, S 2 ) G 5 = 
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Assume that set T is replaced by a superset that is comprised of a finite union of 
rectangles; this replacement does not violate the bound on the truncation error formulated 
in Section [3l Consider a uniform partition (using squared partition sets of size <5) for 
the set T. We employ a bilinear interpolation within each partition set Ai = x 

[ai 2 ,bi 2 ] with basis {(pii{s) = l, 4 >i 2 is) = si,4>i3is) = S2,</>i4('S) = G (or their 

equivalent Lagrange polynomials m)- Assume that the density function ts(-|s) is partially 
differentiable and define the following bounds on its derivatives 


9^ts(5 5 ) 


d'^tfs\s) 

94 

^ yv'i2, 

dsps-k 


<Ml, A: G Na, Vs,s G T. 

The operator IIt in (I5.6p . constructed with bilinear interpolation within each partition set, 
satisfies the inequality 


llnx {ts{-\s)) - ts(-|s)||oo < 


16 


(yAi\ + Ada) + 


,53 

8^2 


(-^3 + -^ 3 ) ) 


Vs G T, 


1 /f? 

where 5 = maxj(5i, with 5i = \{hii — au)^ + {bi 2 — 012 )^] , i G Nn- We implement Algo¬ 

rithm [T] for two-dimensional processes using bilinear interpolation: the approximation error 
is upper-bounded by the quantity 


ll/it - V’t Iloo < e\ = £K{t,Mf), t G Njv, 

with the constant defined in (|5.3p and computed for this particular choice of basis 
functions and points. It can be proved that for bilinear interpolation basis functions, the 
constant of (15..dj) is upper bounded by Mf of Assumption 12.31 and thus can be replaced 
by this quantity in the error computation. 


6.4. Trilinear Interpolation for Three-Dimensional Systems. We now apply the re¬ 
sults of Section [5] to a general three-dimensional Markov process, where s = (si, 52 , 53 ) G 
5 = Again we replace the set T by a superset that is comprised of a finite union of 
boxes, without violating the bound on the truncation error formulated in Section [3l Con¬ 
sider a uniform partition (using cubic sets of size 6) for the set T. We employ a trilinear 
interpolation within each partition set with basis functions 


= {1, 'Sl, 52, 53, 5 i 52, 5253, 5353, 535253}. 


Assume that the density function fa(-15) is partially differentiable and define the following 
bounds on its derivatives 


d‘^ts{s\s) 

< A/f* 

d^tfs\s) 

<r' 

d'^ts{s\s) 

dsf 


dsjsj 

— ■'^'■3 5 

553S2S3 


< Ads, fjjGNs, 

We implement Algorithm [T] for three-dimensional processes using trilinear interpolation in 
the operator IIx ()5.6I) . The approximation error is then upper bounded by the quantity 


WlJ-t -'iptWoo < t G Nat, 

with the constant 

^ (Ml + Ml + Ml) + ^ {Ml^ + Mf + Mf + Mf + Mf + Mf + 3M3) • 


Similar to the bilinear interpolation case, the constant can be replaced by Mf of As¬ 
sumption [2]3] in the error computation. 
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7. Application of the Formal Approximation Procedure 
TO the Probabilistic Invariance Problem 

The problem of probabilistic invariance (or, equivalently, safety) for general Markov pro¬ 
cesses has been theoretically characterised in [3] and further investigated computationally 
in [SmolEIKIl]. With reference to a discrete-time Markov process over a continuous 
state space 5, and to a safe set A € ^(5), the goal is to quantify the probability 

(A) = P{s(t) G A, for all t € Z 7 v|'S( 0 ) = s}. 

More generally, it is of interest to quantify the probability (A), where the initial condition 
of the process s(0) is a random variable characterised by the density function 7ro(-). In 
Section rrn we present a forward computation of probabilistic invariance by application of 
the approximation procedure above, then review results on backward computation laiini 
dniii] in Section EH We conclude in Section o with a comparison of the two approaches. 


7.1. Forward Computation of Probabilistic Invariance. The technique for approxi¬ 
mating the density function of a process in time can be easily employed for the approximate 
computation of probabilistic invariance. Define functions Wt : S —)■ M-^, characterised as 

lW+i(s) = l^(s) y Wi(s)4(s|s)ds, VFo(s) = l^(s)7ro(s), Vs G 5. (7.1) 

Then the solution of the problem is obtained as W]y(s)ds. A comparison of 

the recursions in (j7.1|) and in (j3.4p reveals how probabilistic invariance can be computed 
as a special case of the general approximation procedure in this work. In applying the 
procedure, the only difference consists in replacing set T by the safe set A, and in restricting 
Assumption 12.31 to hold over the safe set - the solution over the complement of this set 
is trivially known, as such the error related to the truncation of the state space can be 
disregarded. The procedure consists in partitioning the safe set, in constructing the Markov 
chain as per m, and in computing as an approximation of Wt{-) based on (14.2p . 
The error of this approximation is \\Wt — il^tWoo < Et, which results in the following: 


I Ms)ds 


IA 


< EnC{A) = k{N,M^)X<^ 5C{A) = Ef. 


Note that the sub-density functions satisfy the inequalities 


1 > 


[ Wo{s)ds > [ Wi{s)d, 

Js Js 


s > ■ ■ ■ > / WN{s)ds > 0. 


7.2. Backward Computation of Probabilistic Invariance. The contributions in [21 
ttomiKii] have characterised specifications in PCTL with an alternative formulation based 
on backward recursions. In particular, the computation of probabilistic invariance can be 
obtained via the value functions 14 : 5 —)■ [0,1], which are characterised as 

f4(s) = Ia('S) j Vt+i{s)ts{s\s)ds, VAr(s) = 1 a(s), Vs G 5. (7.2) 

The desired probabilistic invariance is expressed as p^^{A) = Vo(s)'7ro(s)ds. The value 
functions always map the state space to the interval [0,1] and they are non-increasing, 
namely Vt(s) < V)+i(s) for any fixed s G 5. The contributions in [21 ITOl [HI |T2] discuss 
efficient algorithms for the approximate computation of the quantity p^^(A), relying on 
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different assumptions on the model under study. The easiest and most straightforward 
procedure is based on the following assumption [2]. 

Assumption 7.1. The conditional density function of the process is globally Lipschitz con¬ 
tinuous with respect to the conditional state within the safe set. Namely, there exists a finite 
constant At,, such that 

|t(s|s) — t(s|sOl ^ s', s G A. 

A finite constant is introduced as = sup 5 g_ 4 ts(s|s)(is < 1. 

The procedure introduces a partition of the safe set A = and extends it to 

S = with An-i-i = 5\AI. Then it selects arbitrary representative points Si G Ai 

and constructs a finite-state Markov chain over the finite state space {si, S 2 , • • •, Sn-i-i}, 
endowed with transition probabilities 

P{si,Sj) = ts{s\si)ds, P{Sn+l,Sj) = S(^ri+i)j, (7.3) 

for all i G N„,J G The error of such an approximation is |llj : 

E^, = K{N,Mi,)X^,5C{A), 

where 6 is the max partitions diameter, and C{A) is the Lebesgue measure of set A. 

7.3. Comparison of the Two Approaches. We first compare the two constructed Markov 
chains. The Markov chain .^p obtained with the abstraction from the forward approach is a 
special case of the Markov chain from the backward approach: in the latter case in fact 
the representative points can be selected intelligently to determine the average probability 
of jumping from one partition set to another. More specifically, the quantities (14.ip are a 
special case of those in (17.311 (based on the mean value theorem for integration). We will 
show that this leads to a less conservative (smaller) error bound for the approximation. 

The forward computation is in general more informative than the backward computa¬ 
tion since it provides not only the solution of the safety problem in time, but also the state 
distribution over the safe set. Further the forward approach may provide some insight to 
the solution of the infinite-horizon safety problem |28ll30j for a given initial distribution. As 
discussed in m, solution of the infinite-horizon safety problem depends on the existence 
of absorbing subsets of the safe set. The outcome of the forward approach can provide 
evidence on the non existence of such subsets. Finally, the forward approach presented 
in Sections [3]l6] for approximating density functions can be used to approximate the value 
functions in the recursion dEU) over unbounded safe sets since we do not require the state 
space (thus also the safe set) to be bounded, while boundedness of the safe set is required 
in all the results in the literature that are based on backward computations. 

Next, we compare errors and related assumptions. The error computations rely on 
two different assumptions; the Lipschitz continuity of the conditional density function with 
respect to the current state or to the next state, respectively. Further, the constants Mf and 
Mb are generally different and play an important role in the form of the error. Mb represents 
the maximum probability of remaining within a given set, while Mf is an indication of the 
maximum concentration of the process evolution towards one state, over a single time-step. 
Mb is always less than or equal to one, while Mj could be any finite positive number. 
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Figure 4: Approximate solution of the probabilistic invariance problem (thin black line), 
together with error intervals of forward (blue band) and backward (red band) 
approaches, for a = 1.2 (left) and a = 0.8 (right). 

Example 12.11 (Continued). The constants ApM^ and Ab,M[, for the one dimensional 
dynamical system of Example 12.11 are 

"^1 ^ o ^ ^ Mi, < 1. 

If 0 < a < 1, the system trajectories converge to an equilibrium point (in expected value). 
In this case the model solution has higher chances of ending up in a neighbourhood of the 
equilibrium in time, and the backward recursion provides a better error bound. If a > 1, the 
system trajectories tend to diverge with time. In this case the forward recursion provides a 
much better error bound, compared to the backward recursion. 

For the numerical simulation we select a safety set A = [0,1], a noise level a = 0.1, and 
a time horizon N = 10. The solution of the safety problem for the two cases a = 1.2 and 
o = 0.8 is plotted in Figure [H We have computed constants Xf = 24.20, Mj, = 1 (in both 
cases), while At, = 29.03, Mf = 0.83 for the first case and Af, = 19.36, Mj = 1.25 for the 
second case. We have selected the center of the partition sets (distributed uniformly over 
the set Al) as representative points for the Markov chain In order to compare the two 
approaches, we have assumed the same computational effort (related to the same partition 
size of (5 = 0.7 x 10“^), and have obtained an error Ef = 0.008, Ef, = 0.020 for a = 1.2 and 
Ej = 0.056, Eb = 0.014 for a = 0.8. The simulations show that the forward approach works 
better for a = 1.2, while the backward approach is better suitable for a = 0.8. Note that the 
approximate solutions provided by the two approaches are very close: the difference of the 
transition probabilities computed via the Markov chains Alp Ala are in the order of 10“®, 
and the difference in the approximate solutions (black curve in Figure 0]) is in the order of 
10“®. This has been due to the selection of very fine partition sets that have resulted in 
small abstraction errors. □ 

Remark 7.2. Over deterministic models, |25j compares forward and backward reachability 
analysis and provides insights on their differences: the claim is that for systems with signifi¬ 
cant contraction, forward reachability is more effective than backward reachability because 
of numerical stability issues. On the other hand, for the probabilistic models under study, 
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the result indicates that under Lipschitz continuity of the transition kernel the backward 
approach is more effective in systems with convergence in the state distribution. If we treat 
deterministic systems as special (limiting) instances of stochastic systems, our result is not 
contradicting with [25] since the Lipschitz continuity assumption on the transition kernels 
of probabilistic models does not hold over deterministic ones. D 

Motivated by the previous example, we study how the convergence properties of a 
Markov process are related to the constant Mp 

Theorem 7.3. Assume that the initial density function 7ro(s) is bounded and that the con¬ 
stant Mf is finite and Mj < 1. If the state space is unbounded, the sequence of density 
functions {7rf(s)|t > 0} uniformly exponentially converges to zero. The sequence of proba¬ 
bilities P{s(t) G .4.} and the corresponding solution of the safety problem for any compact 
safe set A exponentially converge to zero. 

Theorem 17.31 indicates that under the invoked assumptions the probability “spreads 
out” over the unbounded state space as time progresses. Moreover, the theorem ensures the 
absence of absorbing sets [281130] , which are indeed known to characterise the solution of 
infinite-horizon properties. Example 17.41 studies the relationship between constant Mj and 
the stability of linear stochastic difference equations. 


Example 7.4. Consider the stochastic linear difference equations 

s{t-\-1) = As{t)-I w{t), s(-), !(;(•) G 

where w{-) are i.i.d. random vectors with known distributions. For such systems = 
l/|detj4|, then the condition My < 1 implies instability of the system in expected value. 
Equivalently, mean-stability of the system implies My > 1. Note that for this class of 
systems Mf > 1 does not generally imply stability, since det4 is only the product of the 
eigenvalues of the system. □ 


The Lipschitz constants Xf and A(, have a different nature, as clarified in Example 17.51 


Example 7.5. Consider the dynamical system 

s{t -h 1) = f{s{t),w{t)), s(-),w(-) G 


where w(-) are i.i.d. with known distribution Suppose that the vector field / : 

X ^ is continuously differentiable and that the matrix ^ is invertible. Then 
the implicit function theorem guarantees the existence and uniqueness of a function g : 

X such that w{t) = g{s{t -|- l),s{t)). The conditional density function of the 

system in this case is [26] : 


ts{s\s) 


det 


dg,. . 


tw{g{s,s)). 


The Lipschitz constants Ap At, are specified by the dependence of function g{s, s) from the 
variables s, s, respectively. As a special case the invertibility of ^ is guaranteed for systems 
with additive process noise, namely /(s, w) = fa{s) + w. Then g{s, s) = s — fa{s), Xf is the 
Lipschitz constant of tu]{-), while At, is the multiplication of the Lipschitz constant of 
and of fa{-). □ 
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8. Conclusions 

This contribution has put forward new algorithms, based on Markov chain abstractions, 
for the efficient computation of approximate solutions of the state distribution function in 
time of Markov processes evolving over continuous state spaces. A higher-order function 
approximation method has also been presented, with a formal derivation of an upper bound 
on the associated error. The approach has been applied to the verification of a particular 
non-nested PCTL formula (expressing probabilistic safety or invariance), and compared 
with an alternative computational approach from the literature. 

The authors plan to integrate the presented procedures within the software tool FAUST^, 
which is developed for the formal abstraction and verification of uncountable-state stochas¬ 
tic processes |16] . The software enables the user to automatically export the hnite-state 
abstracted model to existing probabilistic model checking tools, such as PRISM and MRMC 
[niEB], for further quantitative analysis, verification, and synthesis objectives. 
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Appendix A. Proof of Statements 


Proof of Theorem \3.1\ We prove the theorem inductively. The statement is trivial for t = 0 
based on Assumption l2.2l Suppose it is true for t, then we prove it for t+1. Take s G 5\At_|_i, 
then 

7rt+i(s) = / ts{s\s)TTt{s)ds = / ts{s\s)7i:t{s)ds + / ts{s\s)7i:t{s)ds. 

Js Js\At JAt 

The first integral is upper bounded by £tM^. The domain of the second integral implies 
that (s, s) G At X S. Combining this with the definition (13.21) of A^+i and s ^ Aj+i results 
in (s, s) ^ r. Then 


'At 


ts(s\s)7rt(s)ds < e / '^t{s)ds < e. 


'At 


Then we obtain 7ri+i(s) < etMf + e = et+i, for all s G 5\At+i. □ 

Proof of Theorem \3.3\ The initial density function /tq satishes the following inequality: 

Iko - Molloo = IKo - IaoTToIIoo = l|l5\Ao7ro||oo = sup{7ro(s), s G 5\Ao} < eo- 
Suppose satishes the inequality. We prove that it is also true for ^t+i- Take any s G 5\T, 
5\T C 5\At+i ^ |vrt+i(s) - /it+i(s)| = vrt+i(s) < £t+i. 

Take any s G T, we have 

|7rt+i(s) - |Ut+i(s)| < / ts{s\s)\'Kt{s) - Ht{s)\ds < £t / ts{s\s)ds < £tM^ < £t+l. 

Js Js 

Proof of Lemma \4^ For any t G N and s, s' G S, 

|7rt(s) - 7rt(s')| < / 7rt_i(s) |ts(s|s) - ts(s'|s)| ds < Af||s - s'll / TTt{s)ds = Xf\\s - s 

Js Js 

Proof of Theorem 14-21 We use the triangle inequality as 

11^4 V’tiloo ^ ll^t /I-tlloo T Wht V^llloo ^ S't T Wht V’tlloo- 
Dehne the set of Lebesgue integrable functions V and its subset Vi: 


□ 


□ 


V = {f :S 


>0 


f{s)ds < oo L 


Pi = / G P 


and the operators 


ts{s\s)f{s)ds < oo, for all s G 5 


P(/)(s) = 1 t(s) J ts{s\s)f{s)ds, 
^ [a f{s)ds 

Ua-.V^Vi, na{f){s) = Y,-^WT^^AM- 


2=1 




Then ^,t{s) is formulated as t—times application of the operator TZ to the initial function 

7ro(s), 

^xt+l = T^iht) ht = Vt G N. 
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Define functions ijjt by the recursive equation: 

lAt+i = (7^a7^)(V’t) ^ lAt = (7^a7^)'(V'o) Vt G N, 

initialised with i/jq = TZaino). Notice that the functions ipt, t G'N, are all piecewise constant 
due to their recursive dehnition. We restrict our attention to the set T since the supports 
of both functions are included in T. The goal is to prove that satisfies and 

ll/^t — V'tiloo < Et- We achieve this goal using induction. The initial function ipQ is of the 
form 


V'o(s) = 7^a(7^o)(s) = ^ 


” fAi'^oiv)dv^ . Poii) 


lpt+l{s) = (TZaTZ){'lpt){s) = 


£{A,) 


■j) J Aj 


for a: 

Tl{ijjt)iv)dv 




. 1 ^C{Ai 

2 = 1 2=1 

Suppose the statement is true for t . Then for any s & Aj, j G 

1 




£{A 


■j) JAj JS 


ts{v\u)iptiu)dudv = ^ j I ts{v\u)'tptiu)dudv 

V J / 2=1 '^'^3 ^ 


1 ^ 


Pt{i) 


E{Aj) E(Ai) J^. 


ts(vlu)dudv = 


1 


We have proved (14.21) . The approximation error for t = 0 is computed using (|2.3I) in 
Assumption 12.31 Take any state s G Ai, i G Nn: 


\Pois) - '!/>o(s)| = 


7ro('S) - 


lAiMv)dv 


C{A, 


^lA,\Ms)-Mv)\dv 

- -£(T)-SV« = -Eo. 


For t > 1 we can write 

ll^f+l -1pt+l\\oo = WEint) - ('^a)^)(V’t)l|oo < WEint) -'R,{'lpt)\\oo + ll^(V’t) - {E,a'R,){'lpt 

The error has two terms. The first term is upper bounded as follows 

|'^(A‘t)(s) - )^(V’t)(s)| < 1 t(s) [ ts{s\s)\fj.t{s) -'ipt{s)\ds < Et [ ts{s\s)ds < MfEt. 

Js Js 

Let us focus on the second term: for any arbitrary state s G j G N^, we have 

1 f 


|7^(V’^)(s)-(7^a7^)(V’t)(s)| = 


J ts{s\u)'il>t{u)du- 


E{Aj) 


ts{v\u)^lJt{u)dudv 


ts{s\u) - 


CiAj) 


ts{v\u)dv 


- / . I \ ts { s \ u ) - ts{v\u)\dvdu < I i>t{u)Xf6du < Xf6. 

Js J _4 Js 


The summation of the two upper bounds leads to that in (14.4p . 


□ 
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Proof of Theorem \5.1[ The definition of implies that 'i/’t+i = with 'ipQ = fiQ. 


Then (|5.5I) is true for t = 0 with = 0. Assume that ||//t — V'J’||oo < E^, then for any 
set, 

|/it+i(s) - = |7^T(/it)(s) - ^T7^T('^/’^)(s)l 

< |7^x(/^t)(s) — nT7^T(A‘t)(s)| + |nx7^T(/^t)(s) — IIxT^xC^/'i )(s)|- 

The first term can be upper bounded based on the linearity of the operator fix as 


|^T(i^t)(s) — nx7^x(Mt)('S)l = 


/X 


ts{s\s)iJt{s)ds - / nx(ts(s|s))/xt(s)ds 

Jr 


< [ iJ.tis)\tg{s\s) - nx(ts(s|s))|ds < T*’ / Ht{s)ds < T*’. 
Jr Jr 

On the other hand, the second term is upper bounded as follows: 

I^x7^x(/X^)(s) - ^x7^x(^/'^)(s)l < f |nx(4(s|s))||/^t(s) - V't ('S)|ds 


<E'P 


□ 


|nx(ts(s|s))|ds < E^M^. 

The addition of the two bounds leads to the statement. 

Proof of Theorem \7.S\ Based on the recursion over the density functions, we have that 

vrt+i(s) < sup{7rt(s), s e S} j ts{s\s)ds < Mj sup{7rt(s), s E 5} 

=> < sup{7ro(s), s E 5}Mp Vf E N. 

Then {tt*} converges to zero uniformly exponentially, with a rate Mf. With focus on the 
safety problem we have P{s(u) E A for all u E N*} < P{s(t) E ^}, where 

P{s(f) E ^} = / 7rt(s)(is < £(^) sup 7ro(s)M|. 

J s^S 

Note that if lim P{s(t) E 5} = lim 1 = 1, then the state-space cannot be bounded under 

t—¥00 t—¥00 


the assumptions. 


□ 
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N = {1,2,3,...} 

Nm = {1,2,... ,m} 

’Em = {0,1,2,... ,m} 

S 

B{S) 

T, 

V 

ts{s\s) 

TTo : 5 ^ 

TTi : 5 ^ 

P 

s(t) = [s(0),s(l),... ,s(t)] 

a, b, a, a 
w(-) 

[/3o,7o] C M 
(^a(u) 

rcS^AoCS 

e 

eo 

£t 

Ao 

Af 

Mf 

II * lloo 
Ia(-) 

At 

T = uiloAt 
S(.) 

At = [A,7t] 

T = A 

n 

5 

An+i = S\T 

p = m 

^(n+l)j 

P{-) 

Po = [po(l),---,Po(« + l)] 


Appendix B. List of Symbols 

the set of natural numbers 
the finite set of natural numbers 
the finite set of non-negative numbers 
discrete-time Markov process 
state space of the Markov process 
Borel cj-algebra on the space S 
stochastic kernel 

probability measure on S for one-step transition of the process 

conditional density function of the process 

density function of the initial state 

density function of s{t) 

probability measure on the product space 

bold typeset is employed for vectors 

Markov chain for forward computation of the density function 
parameters of the example 
process noise 

support of the initial density function tto in the example 
Gaussian density function with zero mean and standard devia¬ 
tion a 

sets used in the truncation procedure 

threshold on truncation part of tgi t<i(s|s) < e for all (s,s) £ 

s^\r 

threshold on truncation part of ttq: 7ro{s) < Eq for all s G 5\Ao 
threshold on truncation part of vr*: 7rt{s) < £t for all s G S\kt 
Lipschitz constants of tto{s) 

Lipschitz constants of G(s|s) with respect to s 
upper bound for the quantities J^4(s|s)ds, s G S 
infinity norm 

indicator function of set A 

support set of 7rt{-) obtained via a recursive procedure 
truncated state space 
set-valued map S : S ^ 2^ 

constant equal to t for = 1 and to (1 — — Mf) for 

Mj A 1 

approximate density function after state space truncation 

support set of 7rt{-) in the example 

selected partition 

partition size 

partition diameter 

partition set associated with absorbing state of Markov chain 
transition probability matrix of the Markov chain 
Kronecker delta function, equal to one for j = (n + 1) and zero 
otherwise 

Lebesgue measure of a set 
pmf of a Markov chain at t = 0 
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Pt = \pt{l),---,pt{n + l)] 

'fpt 

Et 

B(5) 

]B(T) 

TZ-f 

^ = {Ms), ■ ■ • 

'I' = span^ 

Ht : B(T) ^ ^ 

4 

E^t 

d 

Ux) 

f = [fisi)]imh 

a = [aj]je% 

Q ~ [4^j{Si)]iJ 
{4‘ij, 3 £ 

{Sij, j £ N/j} 

) j £ 

Qi — ['/’ij 

fUi 
cx^- ■ 

3^ 

h'uv 

puv 

P= [Pii,3)kj 

at = [MM 
Si 

[ai,bi] 

Mh 

[aii,bii] X [ai2,bi2] 

M^,Ml 

pMA) 

P^oiA) 

Wt:S ^ 


pmf of a Markov chain at time t 

approximation of the density function tt^ after abstraction 

abstraction error related to partitioning 

space of bounded and measurable functions on S 

space of bounded and measurable functions on T 

linear operator on B(T) defined as TZr {f){s) = 

fj. t^(sls)f(s)ds, for all s G T 

set of basis functions (/>j G B(T) 

function space generated by 

projection operator 

upper bound for ||nT(4(-|'S)) - is(-|'S)|loo 

upper bound for fj. Iflx (t 5 (s|s))| ds, for all s G T 

approximation of vrt(-) computed via higher-order methods 

error of higher-order approximation 

dimension of the state space 

interpolation operator projecting any function / : P —E to a 
unique function of such that n-p(f) = J2j=i ^j4>j 
/i-dim. column vector containing values of / at the interpola¬ 
tion points 

/i-dim. column vector containing interpolation coefficients 

/i-dimensional interpolation matrix corresponding to IIxi 

set of basis functions for partition set Ai 

set of interpolation points in the partition set Ai 

set of interpolation coefficients in the partition set Ai 

matrix representation of interpolation inside partition set Ai 

/i-dimensional identity matrix 

function / : T —>■ E with domain restricted to set Ai C T 
interpolation coefficients for used in Algorithm [1] 
defined as = TZ-riiptMisuv) and used in Algorithm [T] 
defined as ts{suv\s)4>ijis)ds and used in Algorithm [1] 

matrix with entries P{i,j) = X 4 . used for piecewise 

constant approximation of the density functions in Algorithmic 
a row vector used in Algorithm [2] for values of ipt 
diameter of the partition set Ai, Si = sup{||s — s'||, s, s' G Mj} 
partition set Ai for one-dimensional systems 
an upper bound for the quantity \d^ts{s\s)/ds^\, for all s, s G T 
partition set Ai for two-dimensional systems 
upper bounds on partial derivatives of in 2 -dim. systems 

upper bounds on partial derivatives of in 3-dim. systems 

safety probability over the set A with time horizon N and initial 
state s 

safety probability over the set A with the initial state admitting 
the density function vro 

sub-density functions for forward computation of the safety 
probability 
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Vf.S^ [ 0 , 1 ] 


Ef = K{N,Mf)Xf5C{A) 
= K{N,M^)X^5CiA) 
Afa 


Mb 


P{Si, Sj) 


value functions for backward computation of the safety proba¬ 
bility 

error of forward computation of the safety probability 
error of backward computation of the safety probability 
Lipschitz constant of t{s\s) with respect to s 
upper bound for t 5 (s|s)ds, for all s G ^ 
finite-state Markov chain obtained via the backward abstrac¬ 
tion approach 

transition probabilities of the Markov chain 
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